Variability of organ-biased genes (Danio rerio)

Setup

Parameters

print(params)

## $input.variability.rds.path
## [1] "../../results/03_Variability_Plots"
## 
## $input.tau.rds.path
## [1] "../../results/10_Expression_Specificity"
## 
## $output.rds.path
## [1] "../../results/12_Organ_Bias"
## 
## $seed.run
## [1] 12345
## 
## $species
## [1] "DRE"
## 
## $species.name
## [1] "Danio rerio"
## 
## $tau.cutoff
## [1] 0.3
## 
## $all.nonzero.matrix
## [1] TRUE
## 
## $randomized.mean.expr
## [1] FALSE

Libraries

library(dplyr)

Functions

source("../functions/helper_functions.R")
source("../functions/expression_specificity_plots.R")
source("../functions/organ_bias.R")
source("../functions/organ_bias_plots.R")
#source("../functions/annotate_boxplots.R")

Session Info

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Big Sur 11.7.8
## 
## Matrix products: default
## LAPACK: /Users/cbucao/miniforge3/envs/fish-ev-gene-function/lib/libopenblasp-r0.3.25.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
##  [1] ggtree_3.0.4         tidytree_0.4.6       treeio_1.16.2        org.Dr.eg.db_3.13.0  GOstats_2.58.0      
##  [6] graph_1.70.0         Category_2.58.0      Matrix_1.5-4.1       AnnotationDbi_1.54.1 IRanges_2.26.0      
## [11] S4Vectors_0.30.2     Biobase_2.52.0       BiocGenerics_0.38.0  BimodalIndex_1.1.9   reshape2_1.4.4      
## [16] tibble_3.1.7         tidyr_1.2.1          cowplot_1.1.1        ggpubr_0.4.0         slider_0.2.2        
## [21] viridis_0.6.3        viridisLite_0.4.1    gplots_3.1.3         ggplot2_3.3.6        matrixStats_0.62.0  
## [26] edgeR_3.34.1         limma_3.48.3         colorBlindness_0.1.9 dplyr_1.0.9         
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.1-0       ggsignif_0.6.4         ellipsis_0.3.2         mclust_6.1.1          
##   [5] XVector_0.32.0         fs_1.6.5               aplot_0.2.3            farver_2.1.1          
##   [9] bit64_4.0.5            fansi_1.0.4            splines_4.1.0          cachem_1.0.8          
##  [13] knitr_1.43             jsonlite_1.8.9         broom_1.0.5            annotate_1.70.0       
##  [17] GO.db_3.13.0           cluster_2.1.6          png_0.1-8              BiocManager_1.30.21   
##  [21] compiler_4.1.0         httr_1.4.6             backports_1.4.1        fastmap_1.1.1         
##  [25] lazyeval_0.2.2         cli_3.6.3              htmltools_0.5.5        tools_4.1.0           
##  [29] gtable_0.3.3           glue_1.8.0             GenomeInfoDbData_1.2.7 Rcpp_1.0.10           
##  [33] carData_3.0-5          jquerylib_0.1.4        vctrs_0.4.1            Biostrings_2.60.2     
##  [37] ape_5.8                nlme_3.1-162           xfun_0.39              stringr_1.5.0         
##  [41] lifecycle_1.0.4        renv_0.17.3            gtools_3.9.4           rstatix_0.7.2         
##  [45] XML_3.99-0.9           MASS_7.3-58.3          zlibbioc_1.38.0        scales_1.2.1          
##  [49] RBGL_1.68.0            yaml_2.3.7             memoise_2.0.1          gridExtra_2.3         
##  [53] ggfun_0.1.7            yulab.utils_0.1.8      sass_0.4.6             stringi_1.7.6         
##  [57] RSQLite_2.3.1          genefilter_1.74.1      highr_0.10             caTools_1.18.2        
##  [61] warp_0.2.0             GenomeInfoDb_1.28.4    rlang_1.1.4            pkgconfig_2.0.3       
##  [65] bitops_1.0-7           evaluate_1.0.1         lattice_0.21-8         purrr_0.3.5           
##  [69] patchwork_1.3.0        labeling_0.4.2         bit_4.0.5              tidyselect_1.2.0      
##  [73] AnnotationForge_1.34.1 GSEABase_1.54.0        plyr_1.8.9             magrittr_2.0.3        
##  [77] R6_2.5.1               generics_0.1.3         DBI_1.1.3              pillar_1.7.0          
##  [81] withr_3.0.2            mgcv_1.8-42            survival_3.5-5         KEGGREST_1.32.0       
##  [85] abind_1.4-5            RCurl_1.98-1.12        crayon_1.5.3           car_3.1-2             
##  [89] KernSmooth_2.23-21     utf8_1.2.3             oompaBase_3.2.9        rmarkdown_2.22        
##  [93] locfit_1.5-9.7         grid_4.1.0             isoband_0.2.7          Rgraphviz_2.36.0      
##  [97] blob_1.2.4             digest_0.6.37          xtable_1.8-4           gridGraphics_0.5-1    
## [101] munsell_0.5.0          ggplotify_0.1.2        bslib_0.5.0

Data

Organ expression specificity (tau)

# Load expression specificity computations from 10_Expression_Specificity.Rmd
## Conditions available in all species
if (!params$all.nonzero.matrix) {
  if (!params$randomized.mean.expr) {
    filtered.samples <- readRDS(
      file.path(params$input.tau.rds.path,
                paste(params$species, "filtered.samples.common.rds", sep="_")))     
  } else {
    filtered.samples <- readRDS(
      file.path(params$input.tau.rds.path,
                paste(params$species, "filtered.samples.common", 
                      paste0("rnd", params$seed.run, ".rds"), sep="_"))) 
  }
 
} else {
  if (!params$randomized.mean.expr) {
    filtered.samples <- readRDS(
      file.path(params$input.tau.rds.path,
                paste(params$species, "filtered.samples.common.nonzero.rds", sep="_")))      
  } else {
    filtered.samples <- readRDS(
      file.path(params$input.tau.rds.path,
                paste(params$species, "filtered.samples.common.nonzero", 
                      paste0("rnd", params$seed.run, ".rds"), sep="_")))     
  }
}

filtered.samples$conditions <- split_gonads_df(filtered.samples$conditions)
expr.rank.first <- filtered.samples$expr.rank.first

Expression variability

# Load processed expression variability estimates from 03_Variability_Plots.Rmd
if (!params$all.nonzero.matrix) {
  expr.var <- readRDS(
    file.path(params$input.variability.rds.path,
              paste(params$species, "jack.resid.cv.bind.rds", sep="_"))) 
} else {
  expr.var <- readRDS(
    file.path(params$input.variability.rds.path,
              paste(params$species, "jack.resid.cv.bind.nonzero.rds", sep="_"))) 
}

expr.var <- split_gonads_df(expr.var)

# Retain only conditions included in expression specificity computations
expr.var.common <- expr.var %>%
  dplyr::filter(Tissue %in% filtered.samples$conditions$Tissue)

Main

Expression specificity (tau)

Plot - Tau distribution

All genes with tau

if (!params$all.nonzero.matrix) {
  plot_tau_distribution(filtered.samples$tau,
                        tau.low=params$tau.cutoff, tau.high=params$tau.cutoff,
                        ylim=c(0.0,1.5))  
} else {
  plot_tau_distribution(filtered.samples$tau,
                        tau.low=params$tau.cutoff, tau.high=params$tau.cutoff,
                        ylim=c(0.0,3.6))    
}

With tau and expression variability data

tau.with.var <- 
  filtered.samples$tau[names(filtered.samples$tau) %in% levels(as.factor(expr.var.common$GeneID))]

if (!params$all.nonzero.matrix) {
  plot_tau_distribution(tau.with.var,
                        tau.low=params$tau.cutoff, tau.high=params$tau.cutoff,
                        ylim=c(0.0,1.5))    
} else {
  plot_tau_distribution(tau.with.var,
                        tau.low=params$tau.cutoff, tau.high=params$tau.cutoff,
                        ylim=c(0.0,3.6))    
}

Variability as a function of tau

Reorder organs for plotting

expr.var.rank.common <- expr.var.common %>%
  dplyr::inner_join(expr.rank.first, by="GeneID")

if ("Sex" %in% colnames(filtered.samples$metadata)) {
  levels.order <- 
    c("brain","eye","intestine","liver",
      "heart","muscle","gills","pectoral_fin",
      "ovary","testis")
} else {
  levels.order <-
    c("brain","eye","intestine","liver",
      "heart","muscle","gills","pectoral_fin","gonads")  
}

# Exclude missing conditions from factor levels
levels.exclude <- setdiff(levels.order, levels(as.factor(expr.var.rank.common$Tissue)))

# Reorder levels for plotting
expr.var.rank.common$Tissue <-
  factor(expr.var.rank.common$Tissue, levels=levels.order, exclude=levels.exclude)

Density plots

Variability rank vs. mean expression

density_plot_variability_vs_mean(expr.var.rank.common, title=params$species.name)

Mean expression vs. tau

density_plot_mean_vs_tau(expr.var.rank.common, title=params$species.name)

Residual variation vs. tau

density_plot_resid_var_vs_tau(expr.var.rank.common, title=params$species.name)

Variability rank vs. tau

density_plot_variability_vs_tau(expr.var.rank.common, title=params$species.name)

Boxplots by tau quantiles

Split by tau quantiles

tau.with.var.quantiles <- quantile(tau.with.var)

expr.var.rank.common$Quantile <-
  apply(expr.var.rank.common, MARGIN=1, function(r) {
    if (r[["tau"]] <= tau.with.var.quantiles["25%"]) { "Q1" }
    else if (r[["tau"]] <= tau.with.var.quantiles["50%"]) { "Q2" }
    else if (r[["tau"]] <= tau.with.var.quantiles["75%"]) { "Q3" }
    else { "Q4" }
  })

print(tau.with.var.quantiles)  

##         0%        25%        50%        75%       100% 
## 0.03910663 0.18645619 0.26382625 0.37485886 0.80821801

Mean expression

boxplot_mean_by_tau_quantiles(expr.var.rank.common, title=species.name)

Residual variation

boxplot_resid_var_by_tau_quantiles(expr.var.rank.common, title=species.name)

Variability rank

boxplot_variability_by_tau_quantiles(expr.var.rank.common, title=species.name)

Variability of organ-biased genes

Classify broadly expressed and organ-biased genes (with subsets)

Tissue <- as.factor(expr.var.rank.common$Tissue)

# Per organ, classify genes by organ bias
if ("Sex" %in% colnames(expr.var.rank.common)) {
  expr.var.bias.by.condition <- vector(mode="list", length=2)
  names(expr.var.bias.by.condition) <- c("F","M")
  
  for (s in c("F","M")) {
    expr.var.bias.by.condition[[s]] <- vector(mode="list", length=length(levels(Tissue)))
    names(expr.var.bias.by.condition[[s]]) <- levels(Tissue)
    
    for (t in levels(Tissue)) {
      expr.var.bias.by.condition[[s]][[t]] <- expr.var.rank.common %>%
        dplyr::filter(Sex==s) %>%
        dplyr::filter(Tissue!=t) %>%
        classify_genes_by_organ_bias(tau.cutoff=params$tau.cutoff)
      
      # Reorder factor levels
      expr.var.bias.by.condition[[s]][[t]]$TissueBias <-
        factor(expr.var.bias.by.condition[[s]][[t]]$TissueBias, levels=c("broad","focal",t,"other"))
    }
  }
} else {
  expr.var.bias.by.condition <- vector(mode="list", length=length(levels(Tissue)))
  names(expr.var.bias.by.condition) <- levels(Tissue)  
  
  for (t in levels(Tissue)) {
    expr.var.bias.by.condition[[t]] <- expr.var.rank.common %>%
      dplyr::filter(Tissue!=t) %>%
      classify_genes_by_organ_bias(tau.cutoff=params$tau.cutoff)
  
  # Reorder factor levels
    expr.var.bias.by.condition[[t]]$TissueBias <-
      factor(expr.var.bias.by.condition[[t]]$TissueBias, levels=c("broad","focal",t,"other"))    
  }
}

Boxplots by organ bias

Mean expression

if ("Sex" %in% colnames(expr.var.rank.common)) {
  for (s in c("F","M")) {
    for (t in levels(Tissue)) {
      print(boxplot_mean_by_organ_bias(expr.var.bias.by.condition[[s]][[t]], 
                                       species=params$species, tissue=t, sex=s))
    }
  }
} else {
  for (t in levels(Tissue)) { 
    print(boxplot_mean_by_organ_bias(expr.var.bias.by.condition[[t]],
                                     species=params$species, tissue=t))    
    }
}

Residual variation

if ("Sex" %in% colnames(expr.var.rank.common)) {
  for (s in c("F","M")) {
    for (t in levels(Tissue)) {
      print(boxplot_resid_var_by_organ_bias(expr.var.bias.by.condition[[s]][[t]], 
                                            species=params$species, tissue=t, sex=s))
    }
  }
} else {
  for (t in levels(Tissue)) {
    print(boxplot_resid_var_by_organ_bias(expr.var.bias.by.condition[[t]],
                                          species=params$species, tissue=t))    
  }
}

Variability rank

if ("Sex" %in% colnames(expr.var.rank.common)) {
  for (s in c("F","M")) {
    for (t in levels(Tissue)) {
      print(boxplot_variability_by_organ_bias(expr.var.bias.by.condition[[s]][[t]], 
                                              species=params$species, tissue=t, sex=s))
    }
  }
} else {
  for (t in levels(Tissue)) {
    print(boxplot_variability_by_organ_bias(expr.var.bias.by.condition[[t]],
                                          species=params$species, tissue=t))    
  }

}

Supplementary

Brain-biased genes

bias <- "brain"

if ("Sex" %in% colnames(expr.var.rank.common)) {
  expr.var.brain.bias <- 
    dplyr::bind_rows(expr.var.bias.by.condition[["F"]][[bias]], expr.var.bias.by.condition[["M"]][[bias]])
  
} else {
  expr.var.brain.bias <- expr.var.bias.by.condition[[bias]]
}

# Same, just different colors
if (params$randomized.mean.expr) {
  boxplot_variability_example(expr.var.brain.bias, tissue=bias, species=params$species.name, color="#FFFAFA")
} else {
  boxplot_variability_example(expr.var.brain.bias, tissue=bias, species=params$species.name, color="#67005D")
  boxplot_variability_example(expr.var.brain.bias, tissue=bias, species=params$species.name, color="#00BFC4")    
}

Liver-biased genes

bias <- "liver"

if ("Sex" %in% colnames(expr.var.rank.common)) {
  expr.var.liver.bias <- 
    dplyr::bind_rows(expr.var.bias.by.condition[["F"]][[bias]], expr.var.bias.by.condition[["M"]][[bias]])
  
} else {
  expr.var.liver.bias <- expr.var.bias.by.condition[[bias]]
}

# Same, just different colors
if (params$randomized.mean.expr) {
  boxplot_variability_example(expr.var.liver.bias, tissue=bias, species=params$species.name, color="#FFFAFA")  
} else {
  boxplot_variability_example(expr.var.liver.bias, tissue=bias, species=params$species.name, color="#FF8A00")
  boxplot_variability_example(expr.var.liver.bias, tissue=bias, species=params$species.name, color="#F8766D")  
}

Ovary- or gonads-biased genes

if ("Sex" %in% colnames(expr.var.rank.common)) {
  bias <- "ovary"
  expr.var.ovary.bias <- 
    dplyr::bind_rows(expr.var.bias.by.condition[["F"]][[bias]], expr.var.bias.by.condition[["M"]][[bias]])
  
  if (params$randomized.mean.expr) {
    boxplot_variability_example(expr.var.ovary.bias, tissue=bias, species=params$species.name, color="#FFFAFA")
  } else {
    boxplot_variability_example(expr.var.ovary.bias, tissue=bias, species=params$species.name, color="#67005D")
    boxplot_variability_example(expr.var.ovary.bias, tissue=bias, species=params$species.name, color="#00BFC4")    
  }
  
} else {
  bias <- "gonads"
  expr.var.gonads.bias <- expr.var.bias.by.condition[[bias]]
  
  if (params$randomized.mean.expr) {
    boxplot_variability_example(expr.var.gonads.bias, tissue=bias, species=params$species.name, color="#FFFAFA")
  } else {
    boxplot_variability_example(expr.var.gonads.bias, tissue=bias, species=params$species.name, color="#67005D")
    boxplot_variability_example(expr.var.gonads.bias, tissue=bias, species=params$species.name, color="#00BFC4")    
  }
}

Examples

Brain-biased genes expressed in the pectoral fin

bias <- "brain"
focal <- "pectoral_fin"
  
if ("Sex" %in% colnames(expr.var.rank.common)) {
  sex <- "F"
  brain.biased.in.pecfin <- 
    expr.var.bias.by.condition[[sex]][[bias]] %>%
    dplyr::filter(Tissue==focal)

} else {
  brain.biased.in.pecfin <- 
    expr.var.bias.by.condition[[bias]] %>%
    dplyr::filter(Tissue==focal)
  
}

if (params$randomized.mean.expr) {
  print(boxplot_variability_example(
    brain.biased.in.pecfin, species=params$species.name, tissue=bias, ylab="Variability rank", color="#FFFAFA"))    
} else {
  print(boxplot_variability_example(
    brain.biased.in.pecfin, species=params$species.name, tissue=bias, ylab="Variability rank", color="#67005D"))
    
  print(boxplot_variability_example(
    brain.biased.in.pecfin, species=params$species.name, tissue=bias, ylab="Variability rank", color="#00BFC4"))    
}

Liver-biased genes expressed in the intestine

bias <- "liver"
focal <- "intestine"
  
if ("Sex" %in% colnames(expr.var.rank.common)) {
  sex <- "F"
  liver.biased.in.intestine <- 
    expr.var.bias.by.condition[[sex]][[bias]] %>%
    dplyr::filter(Tissue==focal)

} else {
  liver.biased.in.intestine <- 
    expr.var.bias.by.condition[[bias]] %>%
    dplyr::filter(Tissue==focal)
    
}

if (params$randomized.mean.expr) {
  print(boxplot_variability_example(
    liver.biased.in.intestine, species=params$species.name, tissue=bias, ylab="Variability rank", color="#FFFAFA"))    
} else {
  print(boxplot_variability_example(
    liver.biased.in.intestine, species=params$species.name, tissue=bias, ylab="Variability rank", color="#FF8A00"))
  
  print(boxplot_variability_example(
    liver.biased.in.intestine, species=params$species.name, tissue=bias, ylab="Variability rank", color="#F8766D"))    
}

Ovary-biased (or gonads-biased) genes expressed in the eye

focal <- "eye"
  
if ("Sex" %in% colnames(expr.var.rank.common)) {
  bias <- "ovary"
  sex <- "F"
  ovary.biased.in.eye <- 
    expr.var.bias.by.condition[[sex]][[bias]] %>%
    dplyr::filter(Tissue==focal)
  
  if (params$randomized.mean.expr) {
    print(boxplot_variability_example(
      ovary.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#FFFAFA"))      
  } else {
    print(boxplot_variability_example(
      ovary.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#67005D"))
  
    print(boxplot_variability_example(
      ovary.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#00BFC4"))     
  }
  
} else {
  bias <- "gonads"
  gonads.biased.in.eye <- 
    expr.var.bias.by.condition[[bias]] %>%
    dplyr::filter(Tissue==focal)

  if (params$randomized.mean.expr) {
    print(boxplot_variability_example(
      gonads.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#FFFAFA"))  
  } else {
    print(boxplot_variability_example(
      gonads.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#67005D"))
  
    print(boxplot_variability_example(
      gonads.biased.in.eye, species=params$species.name, tissue=bias, ylab="Variability rank", color="#00BFC4"))  
  }
}

Multiple pairwise comparisons

Description

Variability of the set of x-biased genes expressed in organ i (not x) vs. variability of other organ-biased genes in organ i (excluding i-biased genes)

e.g., variability of brain-biased genes expressed in the ovary vs. variability of other organ-biased genes expressed in the ovary (excluding ovary-biased genes)

Exclude broadly expressed and focal-biased genes

if ("Sex" %in% colnames(expr.var.rank.common)) {
  expr.var.bias.comparisons <-
    lapply(c("F"="F","M"="M"), function(s) {
      dplyr::bind_rows(expr.var.bias.by.condition[[s]], .id="Bias") %>%
        dplyr::filter(!(TissueBias %in% c("broad","focal")))
      }) %>%
    dplyr::bind_rows()  
} else {
  expr.var.bias.comparisons <- expr.var.bias.by.condition %>%
    dplyr::bind_rows(.id="Bias") %>%
    dplyr::filter(!(TissueBias %in% c("broad","focal")))
}

Mean expression

Effect size

comparisons.mean.expr <- list()

comparisons.mean.expr$effect.size <- 
  compute_glass_delta_effect_size(expr.var.bias.comparisons, value="mean")

## `summarise()` has grouped output by 'Bias', 'Tissue', 'Sex'. You can override using the `.groups` argument.

Pairwise Wilcoxon test

# Complete
comparisons.mean.expr$wilcox.complete <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons,
    effect.size.df=comparisons.mean.expr$effect.size,
    value="mean", apply.pvalue.cutoff=FALSE,
    order.levels=levels(Tissue)
  )

# Keep significant observations only
comparisons.mean.expr$wilcox.signif <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons, 
    effect.size.df=comparisons.mean.expr$effect.size,
    value="mean", apply.pvalue.cutoff=TRUE, adj.p=0.05, 
    order.levels=levels(Tissue))

Plot - Heatmap

heatmap_effect_size_pairwise_comparisons_blue2red(
  comparisons.mean.expr$wilcox.complete, species=params$species.name, value="mean", order.levels=levels.order)

heatmap_effect_size_pairwise_comparisons(
  comparisons.mean.expr$wilcox.signif, species=params$species.name, value="mean", order.levels=levels.order)

Log-coefficient of variation

Effect size

comparisons.log2cv <- list()

comparisons.log2cv$effect.size <- 
  compute_glass_delta_effect_size(expr.var.bias.comparisons, value="log2cv")

## `summarise()` has grouped output by 'Bias', 'Tissue', 'Sex'. You can override using the `.groups` argument.

Pairwise Wilcoxon test

# Complete
comparisons.log2cv$wilcox.complete <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons,
    effect.size.df=comparisons.log2cv$effect.size,
    value="log2cv", apply.pvalue.cutoff=FALSE,
    order.levels=levels(Tissue)
  )

# Keep significant observations only
comparisons.log2cv$wilcox.signif <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons, 
    effect.size.df=comparisons.log2cv$effect.size,
    value="log2cv", apply.pvalue.cutoff=TRUE, adj.p=0.05, 
    order.levels=levels(Tissue))

Plot - Heatmap

heatmap_effect_size_pairwise_comparisons_blue2red(
  comparisons.log2cv$wilcox.complete, species=params$species.name, value="log2cv", order.levels=levels.order)

heatmap_effect_size_pairwise_comparisons(
  comparisons.log2cv$wilcox.signif, species=params$species.name, value="log2cv", order.levels=levels.order)

Residual variation

Effect size

comparisons.resid.var <- list()

comparisons.resid.var$effect.size <- 
  compute_glass_delta_effect_size(expr.var.bias.comparisons, value="residual_variation")

## `summarise()` has grouped output by 'Bias', 'Tissue', 'Sex'. You can override using the `.groups` argument.

Pairwise Wilcoxon test

# Complete
comparisons.resid.var$wilcox.complete <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons,
    effect.size.df=comparisons.resid.var$effect.size,
    value="residual_variation", apply.pvalue.cutoff=FALSE,
    order.levels=levels(Tissue)
  )

# Keep significant observations only
comparisons.resid.var$wilcox.signif <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons, 
    effect.size.df=comparisons.resid.var$effect.size,
    value="residual_variation", apply.pvalue.cutoff=TRUE, adj.p=0.05, 
    order.levels=levels(Tissue))

Plot - Heatmap

heatmap_effect_size_pairwise_comparisons_blue2red(
  comparisons.resid.var$wilcox.complete, species=params$species.name, value="residual_variation", order.levels=levels.order)

heatmap_effect_size_pairwise_comparisons(
  comparisons.resid.var$wilcox.signif, species=params$species.name, value="residual_variation", order.levels=levels.order)

Variability rank

Effect size

comparisons.expr.var <- list()

comparisons.expr.var$effect.size <- 
  compute_glass_delta_effect_size(expr.var.bias.comparisons, value="variability")

## `summarise()` has grouped output by 'Bias', 'Tissue', 'Sex'. You can override using the `.groups` argument.

Pairwise Wilcoxon test

# Complete
comparisons.expr.var$wilcox.complete <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons,
    effect.size.df=comparisons.expr.var$effect.size,
    value="variability", apply.pvalue.cutoff=FALSE,
    order.levels=levels(Tissue)
  )

# Keep significant observations only
comparisons.expr.var$wilcox.signif <- 
  compute_pairwise_comparisons(
    comparisons.df=expr.var.bias.comparisons, 
    effect.size.df=comparisons.expr.var$effect.size,
    value="variability", apply.pvalue.cutoff=TRUE, adj.p=0.05, 
    order.levels=levels(Tissue))

Plot - Heatmap

heatmap_effect_size_pairwise_comparisons_blue2red(
  comparisons.expr.var$wilcox.complete, species=params$species.name, value="variability", order.levels=levels.order)

heatmap_effect_size_pairwise_comparisons(
  comparisons.expr.var$wilcox.signif, species=params$species.name, value="variability", order.levels=levels.order)

Save

Observed

if (!dir.exists(params$output.rds.path)) { dir.create(params$output.rds.path, showWarnings=FALSE) }

if (!params$all.nonzero.matrix) {
  saveRDS(expr.var.rank.common,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.rank.common.rds", sep="_")))
  
  saveRDS(expr.var.bias.by.condition,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.by.condition.rds", sep="_")))
  
  saveRDS(expr.var.brain.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.brain.bias.rds", sep="_")))
  
  saveRDS(expr.var.liver.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.liver.bias.rds", sep="_")))
  
  saveRDS(expr.var.bias.comparisons,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.comparisons.rds", sep="_")))
  
  saveRDS(comparisons.mean.expr,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.mean.expr.rds", sep="_")))
  
  saveRDS(comparisons.log2cv,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.log2cv.rds", sep="_")))
  
  saveRDS(comparisons.resid.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.resid.var.rds", sep="_")))

  saveRDS(comparisons.expr.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.expr.var.rds", sep="_")))
  
} else {
  saveRDS(expr.var.rank.common,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.rank.common.nonzero.rds", sep="_")))
  
  saveRDS(expr.var.bias.by.condition,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.by.condition.nonzero.rds", sep="_")))
  
  saveRDS(expr.var.brain.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.brain.bias.nonzero.rds", sep="_")))
  
  saveRDS(expr.var.liver.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.liver.bias.nonzero.rds", sep="_")))
  
  saveRDS(expr.var.bias.comparisons,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.comparisons.nonzero.rds", sep="_")))
  
  saveRDS(comparisons.mean.expr,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.mean.expr.nonzero.rds", sep="_")))
  
  saveRDS(comparisons.log2cv,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.log2cv.nonzero.rds", sep="_")))
  
  saveRDS(comparisons.resid.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.resid.var.nonzero.rds", sep="_")))

  saveRDS(comparisons.expr.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.expr.var.nonzero.rds", sep="_")))
}

Randomized mean expression

if (!dir.exists(params$output.rds.path)) { dir.create(params$output.rds.path, showWarnings=FALSE) }

if (!params$all.nonzero.matrix) {
  saveRDS(expr.var.rank.common,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.rank.common", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))

  saveRDS(expr.var.bias.by.condition,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.by.condition", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.brain.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.brain.bias", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.liver.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.liver.bias", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.bias.comparisons,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.comparisons", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.mean.expr,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.mean.expr", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.log2cv,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.log2cv", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.resid.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.resid.var", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))

  saveRDS(comparisons.expr.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.expr.var", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
} else {
  saveRDS(expr.var.rank.common,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.rank.common.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))

  saveRDS(expr.var.bias.by.condition,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.by.condition.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.brain.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.brain.bias.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.liver.bias,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.liver.bias.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(expr.var.bias.comparisons,
          file.path(params$output.rds.path, 
                    paste(params$species, "expr.var.bias.comparisons.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.mean.expr,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.mean.expr.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.log2cv,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.log2cv.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
  
  saveRDS(comparisons.resid.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.resid.var.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))

  saveRDS(comparisons.expr.var,
          file.path(params$output.rds.path, 
                    paste(params$species, "comparisons.expr.var.nonzero", 
                          paste0("rnd", params$seed.run, ".rds"), sep="_")))
}